clear
format long



%--------------------------------------------------------------------------
%----------------------NOTES-----------------------------------------------
%--------------------------------------------------------------------------

%NAME: CRTBPsingleOrbit.m

%Dependencies:
%This program depends on 'CRTBP.m', 'BackwardCRTBP.m', 'magnitudeVelocity.m',
%'librationPoints.m', 'jacobiConst.m', 'interpCrossingData.m', and 
%'CRTBPpoincareNewton.m'.  

%All of these are of these are called by the present program except
%'BackwardCRTBP.m' which is called by 'CRTBPpoincareNewton.m'.

%USE: This is a version of 'CRTBPpoincare.m' which is desinged for looking at 
%specific orbits read off the Poincare map generated by 'CRTBPpoincare.m'
%or one of it's decendants.

%Below are several USER DEFINED CONTROL PARAMETERS.  The user should
%set these (or leave them set to their default values) and save the
%program.  Then this program is typically run from the command prompt.

%OUTPUT:
%The output is a graph of the trajectory and a graph of it's intersection 
%with the Poincare section.

%PURPOSE: After isolating a point on the Poincare section of interest, the
%purpose of this program is to recompute the trajectory of just that point
%and plot that information various ways.




%--------------------------------------------------------------------------
%--------------------USER DEFINED CONTROL PARAMETERS-----------------------
%--------------------------------------------------------------------------

%The parameter 'G' is available in case non-normalized units are
%desired.  When this is not the case just leave 'G' set to one.
        G=1;

%The value of 'mu' depends on the mass ratio fo the primaries.
%If mu is zero the CRTBP reduces to the Kepler problem in non-inertial
%(rotating) coordinates.  The coordinates are normailized so that mu is
%always less than or equal to 1/2, with 1/2 the case of equal masses.
%Other important values are
%
%mu=3.054248395e-06      (for the sun/earth system)
%mu=0.012277471          (for the earth/moon system)
%
%The value of mu must be set by the user (default is earth/moon).
    mu=0.012277471

%Once 'G' and 'mu' are defined we can compute the location of the libration
%points and their Jacobi Integrals.  These will be handy when we set the
%value of the Jacobi Integral for the simulation below.

    %AXULARY CONSTANTS:
    %Compute the location of the five libration points as points in three
    %space.  These can be used to set the value of the Jacobi Constant.
    [L1, L2, L3, L4, L5]=librationPoints(mu);

    %Output location to the screen;
    outL1=L1;
    outL2=L2;
    outL3=L3;
    outL4=L4;
    outL5=L5;
    %Now compute the Jacobi Energy at each collinear libration point

        %The velocity at any equilibrium is zero.
        v_equil=[0; 0; 0];

    %Compute the energy at each libration point
    CL3=jacobiConst(L3, v_equil, mu);
    CL1=jacobiConst(L1, v_equil, mu);
    CL2=jacobiConst(L2, v_equil, mu);
    CL4=jacobiConst(L4, v_equil, mu);
    %CL5=CL4 so no extra computation is needed

%The value of 'C' is the value of the Jacobi Integral that you want to
%compute a Poincare section for.  For values less than CL1 but greater than
%CL2 the 'neck' at L1 is open but the 'neck' at L2 is closed.  Then we can
%have transport between the earth and the moon, but all orbits (in the
%Hill's Region of the Earth and Moon) are bounded.  Then, while the user
%can set 'C' to anything they like, our default value of 'C' is:

    C=CL1-0.5*(CL1-CL2)    %0.5 here is arbitrary.  Any coefficent between
                           %zero and one (inclusive) will have the property
                           %described above.

    %A secondary default value is provided below.  At this energy the neck at
    %L1 is not open.  Hence for initial conditions begining inside the Hills
    %region of the primary, the time scale problems are avoided

    %C=1.001*CL1;   %(Comment out this line if using primary default)

%The variable 'itterates' decides how many crossings of the xz-plane will
%be computed for every initial condition.  This is equivalent to specifying
%how many itterates of the poincare map to compute for every point. Values
%near 1000 are best, but involve long runtimes.
    iterates=1000
    
%These are parameters uniquely determine which orbit will be analized.  
%'PoincareMap1' is the x coordinate from the poincare map of the point 
%you want to examine.  'PoincareMap2' is the xdot or vertical coordinate.
poincare_x=  0.0915;
  
poincare_y=  0.0;


%--------------------------------------------------------------------------
%------------COMPUTATIONS--------------------------------------------------
%--------------------------------------------------------------------------



%Initial Positions
x0=poincare_x;
y0=1e-12;
z0=0.0;

%--------Comment all this out if specifying initial velociety--------------
%Compute the initial velociety from the position and energy and a direction
    
    %Direction of Velociety
    v_x0=poincare_y;
    
   
    %then the magnitude
    V=magnitudeVelocity([x0; y0; z0],C, mu);;
    v_y0=sqrt(V^2 - v_x0^2);
    v_z0=0.0;
    v0=[v_x0; v_y0; v_z0];
    %make it a unit vector
    u_v=v0/norm(v0);
    %then the initial velociety is the unit vector times the magnitude
    xdot0=V*u_v(1,1);
    ydot0=V*u_v(2,1);
    zdot0=V*u_v(3,1);

    
%Initial condition vector
y0=[x0;
   y0;
   z0;
   xdot0;
   ydot0;
   zdot0];
    
    
 %-------------------------------------------------------------------------   
 %All of the variables below do not need to be touched in this program.  
 %They are in some sense irrelevant, but keeping them allows less rewriting 
 %of 'CRTBPpoincare.m'.
 %-------------------------------------------------------------------------   
    
    
%The variable 'k' determins how many initial conditions will be used to
%generate the Poincare map, or equivelently how many points in the domain
%of the mapping to sampled.  For this program k should always be 1.
    k=1




%The variable 'x_begin' sets the left endpoint for the interval of points
%that will be itterated.  Want to stay away from -mu and 1-mu as these are
%singularities.
    x_begin=poincare_x;

%The variable 'x_end' is the right end-point.  The determines how far to
%go to the right of 'x_begin'.
    x_end=x_begin;

%Secondary Controls:
%NOTE: The following parameters do not need to be changed (they were fine
%tuned in the developement of the program.  Changing them should not hurt
%and could even help, but is not necessary).

    %sets the tolerence for the map.  The poincare map looks for crossings
    %of the xy plane.  epsilon determines how close a trajectory point must
    %be to the plane to be considered "on" the plan.  Standard value is
    %ten to the minus six.
        epsilon=1e-6

%The program integrates for a time, then checks for crossings.  It repeats 
%process untill it has found the desired number of crossings.  The variable
%below defines how long to integrate before looking for crossings.  Thsi is
%somewhat arbitrary.

   timeStep=10;
        
        
        
%--------------------------------------------------------------------------



%NOTE: The following parameters are necessary but depend on the above

    %if k=1 no steps are necessare and the following variable makes
    %no sense.
    if k>1
        %determine the step size between successive initial conditons
        xStep=(x_end-x_begin)/(k-1);
    else
        %if k=1 just set step to zero so it is at least defined later
        xStep=0;
    end %end seeing if steps are necessary

%INITIALIZE:
    %This is the number of iterates of the current initial conditions
    %found this far through the while loop
    x_nIterates=0;
    %This is the total number of points in the map so far
    totalIterates=0;
    %Local name of initial conditions:
    initial_n=0;
    %number of times through the while loop
    checkWhile=0;
    lastPath=0;
    pathTime=0;
    %number o times the intorpolator is called
    interpTimes=0;
    ReIntegrate=0;
    checkTotal=0;
    calledNewton=0;
    tf=0;
    
    
%--------------------------------------------------------------------------
%----------------------COMPUTATION OF POINCARE MAP-------------------------
%--------------------------------------------------------------------------

for n=1:k
    %I like to see what loop I'm in at command line, just to see how fast
    %the program is running and where it is.  Comment this out or delete it
    %if you prefer.  Same goes for any 'check' variable
    check_n=n;
    checkTotal=totalIterates;
    %First the position
         
    %Set the initial condition at the begining of the nth time
    %through the for loop (nth initial condition) using the data above:
            initial_n=y0';

    %THE MAIN WHILE LOOP:  This computes the number of times the
    %trajectory with initial condition 'initial_n' crosses the poincare
    %section with positive y velocity. The desired number of such crossings
    %is given by the integer 'iterates'.  The count of how many have been
    %found so far is in the raviable:
              x_nIterates=0;   %This should be set back to zero before
                               %every run through the while loop.
    
    
    %Begining of while loop
    while x_nIterates <iterates
        
               
        %Intigrate this for time 'timeLength'.
           tspan=[0 timeStep];
           %numerical integration
           options=odeset('RelTol',1e-13,'AbsTol',1e-20);     %set tolerences
           [t, trajectory_n]=ode113('CRTBP',tspan,initial_n,options,[],G,mu);

       %how big was the trajectory this time?
       temp=size(trajectory_n);
       numSteps=temp(1,1);    
            
       %keep track of the whole trajectory up to now
       path(lastPath+1:lastPath+numSteps, 1:6)=trajectory_n(1:numSteps,1:6);
       pathTime(lastPath+1:lastPath+numSteps)= timeStep*(checkWhile-1)*ones(1,1:numSteps)+t(1:numSteps);
       
       %incriment the while counter
       checkWhile=checkWhile+1
       %keep track of how much is in path up to now       
       temp2=size(path);
       lastPath=temp2(1,1);
           
           
       %Now check for pairs of points along the trajectory between which
       %the sign of y changes and the sign of ydot is positive
       for i=2:numSteps
           %Check for sign change of y and correct sign of ydot.
           if (sign(trajectory_n(i,2)) ~= sign(trajectory_n(i-1,2))) & (trajectory_n(i,5) > 0) 
                %It could happen that one of these two points is already a
                %crossing to the desired tolerence.  We check for this:
                if (abs(trajectory_n(i-1,2))<epsilon)
                    intersection=trajectory_n(i-1,:);
                elseif (abs(trajectory_n(i,2))<epsilon)
                    intersection=trajectory_n(i,:);
                else
                %If (as is most likely the case) neither of the points 
                %for which the sign change occurs are in tolerence
                %to be considered as crossings and we are farther than 'delta' 
                %from a singularitys, then a Newton method determines the 
                %crossing to the desired accuracy.
                    intersection=CRTBPpoincareNewton(t(i-1),trajectory_n(i-1,:),t(i),trajectory_n(i,:),epsilon,G,mu);
                    calledNewton=calledNewton+1;  %(Print to scren)
                end %end of the intersection finding conditional.

                %Keep track of how many iterates of the current initial
                %condition have been found
                    x_nIterates=x_nIterates+1;  %(print to screne)
                %Keep track of how many total intersections have been found
                    totalIterates=totalIterates+1 %(print to screne)
                %Stor the intersection information for plotting later
                    PoincareMap(totalIterates, 1)=intersection(1);
                    PoincareMap(totalIterates, 2)=intersection(4);
           end %end of pair checking if statement
       end %end of pair checking loop (i=2:numSteps index)

        %make sure variables are ready for next while loop: unless all the
        %iterates have been found we will go fhrought the while loop again.
        %In that case the new initial condition is the final condition of
        %trajectory_n

        initial_n=trajectory_n(numSteps,:);
        tf=t(numSteps)+tf
    end %end of main while loop (x_nIterates <= iterates)

    %The set up for the next time through the 'n' loop is done at the
    %begining of the loop instead of here at the end.
    %So just end and go back.

end  %end of outmost for loop (n=1:k index)

%save trajectoryData1 PoincareMap

figure
hold on

plot(PoincareMap(:,1), PoincareMap(:,2),'.')

figure
hold on
%plot the libration points
plot(L3, 0, 'ko', L1, 0, 'ko', L2, 0,'ko', L4(1,1), L4(2,1),'ko',  L5(1,1), L5(2,1),'ko')

%plot the planets (primaries)
plot(-mu,0, 'r*', 1-mu, 0, 'g*')

%plot the forward trajectory
plot(path(:,1), path(:,2), 'b')



%integrate the problem again to give a whole trajectory
    
%    tspan=[0 tf];
%    %numerical integration
%    options=odeset('RelTol',1e-9,'AbsTol',1e-6);     %set tolerences
%    [t,x]=ode113('CRTBP',tspan,y0,options,[],G,mu);

%figure
%hold on

%plot the libration points
%plot(L3, 0, 'ko', L1, 0, 'ko', L2, 0,'ko', L4(1,1), L4(2,1),'ko',  L5(1,1), L5(2,1),'ko')

%plot the planets (primaries)
%plot(-mu,0, 'r*', 1-mu, 0, 'g*')

%plot the forward trajectory
%plot(x(:,1), x(:,2), 'b')









